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Theoretical and experimental values to date for the resistances of single molecules commonly 
disagree by orders of magnitude. By reformulating the transport problem using boundary conditions 
suitable for correlated many-electron systems, we approach electron transport across molecules from 
a new standpoint. Application of our correlated formalism to benzene-dithiol gives current-voltage 
characteristics close to experimental observations. The method can solve the open system quantum 
many-body problem accurately, treats spin exactly and is valid beyond the linear response regime. 

PACS numbers: 73.63.-b, 73.63.Rt, 71.10.-w, 73.40.Gk 



Substantial interest in recent years has been applied 
to extend the machinery of electronic structure theory, 
largely developed for study of closed or periodic systems, 
to open quantum systems. A primary motivation for this 
effort is the study of electronic transport for nanoscale 
systems, and specifically the study of molecular electron- 
ics P, 0, Recent approaches have introduced den- 
sity functional theory (DFT) into scattering jj, la] or 
Green's function calculations of transport [a, 13 U @- 
The calculated single-electron transmission coefficients 
T are then used to predict a molecular conductance 
using the Landauer-Biittiker formula a = TGo, where 
Go = 2e 2 /h is the conductance quantum |10|]. Notwith- 
standing this effort, the currents flowing through molec- 
ular scale systems are consistently predicted to be sev- 
eral orders of magnitude larger than experimentally re- 
ported values [Hl.llll lj. By formulating scattering bound- 
ary conditions appropriate to many-body wavefunctions, 
we are able to present correlated transport calculations 
which show that an accurate description of electron cor- 
relation is essential for single molecule transport. We find 
current magnitudes close to experiment for the prototyp- 
ical molecular electronics system, benzene-dithiol bonded 
between two gold contacts. 

Central to the study of electronic transport is the con- 
cept of two reservoirs locally in equilibrium, but driven 
out of equilibrium with respect to one another. The cor- 
responding electrochemical potential imbalance is inter- 
preted as the voltage driving a current across a device 
in contact with both reservoirs. Typically, scattering 
boundary conditions are introduced to define the prop- 
erties of the electron reservoirs: one reservoir (the left) 
is associated with the emission of incoming, right mov- 
ing particles with a fixed equilibrium energy distribution, 
but capable of absorbing outgoing, left moving particles 
at any energy; the converse is true for the right reser- 
voir. Implementing scattering boundary conditions for 
a system of fermions within a single particle approxima- 
tion is straightforward. Fermi-Dirac statistics are applied 
to the left and right reservoir incoming wavefunctions 
ipf' R (r) and a voltage is introduced as the difference in 
chemical potential between the left and right reservoir 



distribution functions np(ei ± eV/2) inducing a net cur- 
rent flow. When this scheme is applied with commonly 
used DFT computations, several formal and technical 
problems arise. For example, the exchange-correlation 
functional is current dependent a fact not explored 
to date within DFT transport studies. Even ignoring 
this current dependence, a choice for approximating the 
exchange-correlation functional must be made. Recently, 
it has been pointed out that the currents theoretically 
predicted with DFT methods can vary by over an or- 
der of magnitude simply due to the choice of exchange- 
correlation functional Another drawback to DFT 
transport implementations is the assumption that the 
Kohn-Sham orbitals tpi are physical single particle states 
and give the correct transmission coefficients. Use of the 
Landauer-Biittiker formula beyond the linear response 
regime also introduces a level of approximation which is 
difficult to control. 

In our approach to the electronic transport problem, 
we avoid these issues by working directly with many-body 
wavefunctions ^(risi, . .. ^t^sn) and the exact molecu- 
lar electronic Hamiltonian H. By discovering that a form 
of the single-electron scattering boundary conditions 
can be generalized to many-body wavefunctions, we can 
apply many-body methods to the electronic structure of 
open systems. Here we use the configuration interac- 
tion (CI) formalism which can give accurate solutions to 
quantum many-body problems with correct spin eigen- 
functions and an equal treatment of ground and excited 
electronic states; our approach is equally applicable to 
any other correlated wavefunction method. In CI we ex- 
pand the N-electron wavefunction \& in terms of spin- 
projected Slater determinants 



* = Ci*i + c 2 * 2 



(1) 



with expansion coefficients c M . We use the exact non- 
relativistic electronic Hamiltonian 

h 2 



H 



= l,N 



(2) 
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where U(r) is the attractive potential energy of an elec- 
tron in the Coulomb field of the atomic nuclei. We then 



diagonalize the matrix of the Hamiltonian in the 
basis "J/^ to find the many-body eigenstates and energies. 

Working directly with the N-particle wavefunction W 
removes a direct physical interpretation of single-electron 
wavefunctions ipi(v) and eigenvalues q from the descrip- 
tion of many-electron systems. Therefore, generalization 
of the conventional single-electron scattering boundary 
conditions is not possible, and a formulation of the trans- 
port problem does not exist for correlated many-electron 
systems. We resolve this issue by formulating the scat- 
tering boundary conditions in terms of the Wigner func- 
tion f(q,p). This function has been used as an alterna- 
tive method for applying scattering boundary conditions 
for one-electron problems |l3l |. yielding similar results 
to the conventional boundary conditions. For one elec- 
tron in one dimension, the real function f(q,p) is defined 
by taking the Wigner transform of the density matrix 
p(x,x') = ip* (x')tp(x), given by 

f{l,p)= / drp(q + r/2,q~r/2) exp(-ipr/h), (3) 



and it has the advantage of echoing the properties of a 
classical distribution function (aside for the requirement 
of strictly positive probabilities) with q and p playing 
the role of position and momentum, thus giving a phase- 
space picture of quantum mechanics. For example, ex- 
pectation values of operator functions A(q) of position q 
may be expressed as 

< A(q) >= -L / dpdqA(q)f(q,p) (4) 

27Tft J 

and similarly for functions of momentum p. The Wigner 
function also satisfies a quantum Liouville equation with 
a well-defined classical limit as H — > leading to the 
Boltzmann transport equation. Single-electron scatter- 
ing boundary conditions can then be applied to a domain 
< x < L by requiring that f(0,p > 0) and f(L,p < 0) 
be constrained to their equilibrium values. That is, the 
momentum distribution of the incoming electrons on the 
left and the right is fixed to values characteristic of the 
reservoir. The Schrodinger equation is then solved, sub- 
ject to these boundary conditions, thereby allowing the 
device to find a steady state by varying the values of the 
Wigner function in the interior of the device region and at 
f{0,p < 0) and f(L,p > 0) (the outgoing electrons). We 
note that these boundary conditions break time-reversal 
symmetry as is necessary to generate a current-carrying 
state. 

To find a form of the scattering boundary conditions 
applicable to the many-body problem, we note that 
the same Wigner transform can be applied to the one- 
particle reduced density matrix calculated from a many- 
body wavefunction 

p(r,r') = / dx 2 ... dxjv**(r / s 1) x 2 , • • ■ ,Xjv) 

si 

x* (rsi,x 2 , . . .,xjv) (5) 
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FIG. 1: (top) Geometry of the Aui3-benzene-dithiol-Aui3 
junction, (middle) The Wigner momentum distribution func- 
tion at the left reservoir/device interface at equilibrium (V = 
0). We use Hartree atomic units for p (m = l,c ~ 137). (bot- 
tom) The difference between the finite- and zero-bias Wigner 
functions at the left contact: note the change in scale. The 
constraints have fixed the momentum distribution for mo- 
menta p > 0, while allowing the distribution at p < to vary. 



where x n = r„s„. This generates the Wigner function 
/(q, p) of the correlated wavefunction 'J. To simplify 
matters, we additionally integrate / over planes perpen- 
dicular to the current flow, yielding a one-dimensional 
Wigner function f(q,p) as above, but our approach does 
not depend on this simplification. Our example sim- 
ulation uses the Aui3-benzene-dithiolate-Aui3 contact- 
molecule-contact system to model the transport region, 
shown in fig. 1 (top). Having found the CI ground state 
of the transport region in the absence of voltage bias, 
we calculate the equilibrium values of the Wigner func- 
tion in the left and right contacts fo(qL,p) for p > and 
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fo(qn,p) for p < 0; the left contact values are shown 
in fig. 1 (middle). We then formulate the quantum 
many-body transport problem by making the following 
Ansatz: the appropriate many-body wave function ^> for 
the transport region in the presence of an applied elec- 
tric field £i along the molecular axis minimizes the exact 
energy E =< + e£ J2 n ^"1^ > while simultaneously 
satisfying the open system boundary conditions 

< *|F(q L ,p)|* >= fo(q L ,p) forp > (6) 

and 

< *|F( te ,p)|* >= h{q R ,p) for p < (7) 

with normalisation constraint <3E'| , B>= 1 and constraints 
to enforce current continuity. Here F(q,p) is the hermi- 
tian operator corresponding to the quantity f{q,p)- For- 
mally our Ansatz corresponds to the zero-temperature 
maximum entropy principle on the transport region with 
the boundary conditions specified as system observables 

To numerically implement this scheme, we replace 
the infinite number of constraints implicit in equa- 
tions ©,0 with a finite set of constraints on a grid in 
momentum space. To solve this non-linear constrained 
minimization problem we introduce Lagrange multipliers 
Ai for each constraint, giving us a simpler unconstrained 
minimization problem of the associated Lagrangian func- 
tion 

L=< *|H + e£^z„|4< > -£}Ai < *|C l |* > (8) 

n i 

where the C l are the hermitian operators corresponding 
to the constraint observables. We need to determine the 
correct values of the Lagrange multipliers; these are not 
known a priori, but can be found by an iterative opti- 
mization technique 0] . In figure 1 (bottom) we show the 
difference between a finite bias and equilibrium Wigner 
distribution and verify that the minimization procedure 
respects the boundary constraints. Electron flow results 
from the difference between these distributions, giving a 
net electric current to the right. Once the minimising 
wavefunction VP has been obtained, we find the current 
/ flowing simply by integrating the probability current 
density J(r) over planes perpendicular to the molecular 
axis. 

For our prototype transport calculation, we begin with 
a fully relaxed structure (no strain at zero current) con- 
sisting of the benzene-dithiolate molecule bonded be- 
tween the (111) faces of two gold clusters in C 2 v symme- 
try [l^; see fig. 1 (top). A set of many-body expansion 
functions for this structure was selected usin g th e Monte 
Carlo configuration interaction method [TtI lisj. Exact 
spin coupling is applied and the zero-bias total system 
eigenfunction is a singlet for the ground state. In fig. 



2, we present our computed current- voltage (IV) charac- 
teristics for the benzene-dithiol system. The IV curve has 
been generated by increasing the applied field £ in steps 
and finding the energy-minimising wavefunction at each 
field. The voltage is then calculated by integrating the 
external electric field between the two reservoir/device 
interface planes where the Wigner constraints are ap- 
plied. 
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FIG. 2: Calculated current- voltage characteristics for 
benzene-dithiolate bonded to the (111) faces of two gold clus- 
ters; also shown is the differential conductance 

We find two distinct regions in the resulting IV charac- 
teristics, with an abrupt transition between the two near 
2.5 V. We can trace the origins of this transition by ex- 
amining the polarization behaviour of the isolated A1113- 
benzene-dithiolate-Au 13 system under application of an 
external field. By removing the Wigner constraints we 
decouple the transport region from the gold wires; in fig. 
3 we show the energies of the first five eigenstates as the 
applied field is increased. The energy levels are labelled 
by the irreducible representation of C 2v by which they 
transform at zero field, and by energy ordering within 
this irrep. The system ground state belongs to A\, and 
at finite field mixes with states of B2 symmetry. Sec- 
ond order perturbation theory predicts that the ground 
state and first excited state \1B 2 > will respond slowly 
to the applied field, as they are isolated in energy and 
so have large energy denominators, as we indeed see. In 
contrast, the three states |2Ai>, |3j4i> and \2B 2 > can 
couple strongly, and a state associated with these exci- 
tations polarizes rapidly with the external electric field, 
and crosses the ground state at approximately 2.5 V. Fur- 
thermore, the rigorous definition of correlation involves 
physics beyond that describable by a single Slater deter- 
minant. The eigenvalues of the reduced density matrix 
p(r, r') are the natural orbital occupancies rii\ a necessary 
condition for ^ to originate from some single determinant 
is that ni — 2 for N/2 occupied orbitals and = for 
all others. The CI wavefunction at V = has a typical 
pattern of occupancies slightly changed from that of a 
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FIG. 3: Energies of the lowest five CI eigenstates as a function 
of applied field, for the isolated (no Wigner constraints) Aui3- 
benzene-dithiolate-Aui3 cluster. 



currents orders of magnitude lower than other theoret- 
ical predictions, but of the correct order of magnitude 
with respect to the experimental measurements, the first 
such result to our knowledge. We hypothesize that as 
DFT exchange correlation functionals have been devel- 
oped to reproduce an integrated quantity (the energy), 
they must be modified for accurate treatment of trans- 
port by determining exchange-correlation potentials that 
are locally accurate. We note that for accurate molecular 
electronics simulation, it appears that a highly correlated 
treatment of the electronic structure is needed. 
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single determinant, and this remains the same up to the 
transition. Above the transition, however, we see a sig- 
nificant number of orbitals whose occupancies differ from 
2 or 0. Taken together with the excited state influence 
on the transition voltage, we see that electron correla- 
tion is essential for an accurate description of molecular 
electronic transport. 

We note here that the current magnitude below reso- 
nance is in good agreement with experiment 0, partic- 
ularly when compared to other theoretical predictions, 
and the qualitative behaviour of the IV characteristic is 
consistent with other recent findings of molecular scale 
transport, with non-resonant tunnel regimes followed by 
a rapid increase in current magnitude Our calcu- 

lations indicate an onset resistance of 18 Mf2 and a res- 
onant tunnel resistance of several Mf2, close to the cor- 
responding experimental values of « 22 MSI and w 13 
Mf2 estimated by Reed et al. Previous ab initio DFT 
calculations typically find currents in the range of 20 000 
to 60 000 nAmps at 2 V 0, Q , whereas we find a cur- 
rent of 160 nAmps which compares favourably with the 
experimental measurement of rs 60 nAmps 0. 

In conclusion, we have presented an approach to elec- 
tronic transport that represents a radical departure from 
either scattering theory or Green's function methods. 
The formalism is conceptually simple, and is made pos- 
sible by a generalisation of the scattering boundary con- 
ditions to many-body calculations. Our method yields 
many-body wavefunctions which are exact eigenstates of 
spin and so can also be used for studying spin-dependent 
transport and Kondo effect physics; we can also easily 
add a gate field to the Hamiltonian. For the AU13- 
benzene-dithiolate-Aui3 system, we calculate electronic 
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